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ABSTRACT 



The BLLac object S4 0954+65 is one of the main targets of the Urumqi monitoring program targeting IntraDay Variable (IDV) 
sources. Between August 2005 and December 2009, the source was included in 41 observing sessions, carried out at a frequency of 
4.8 GHz. The time analysis of the collected light curves, performed by applying both a structure function analysis and a specifically 
developed wavelet-based algorithm, discovered an annual cycle in the variability timescales, suggesting that there is a fundamental 
contribution by interstellar scintillation to the IDV pattern of the source. The combined use of the two analysis methods also revealed 
that there was a dramatic change in the variability characteristics of the source between February and March 2008, at the starting time 
of a strong outburst phase. The analysis' results suggest that the flaring state of the source coincides with the appearance of multiple 
timescales in its light curves, indicating that changes in the structure of the relativistically moving emitting region may strongly 
influence the variability observed on IDV timescales. 

Key words, quasars: individual: S4 0954+65; scattering; radio continuum: galaxies; methods: data analysis 



1. Introduction 

Radio variability on timescales from days to weeks was observed 
in AGNs as early as 1971 (see, e.g., Kinman & Conklin [197T] 
Wills [19711 In 1986, Witzel et al. ( fW86l l reported the discovery 
of flux density variations occurring on timescales shorter than 
one day, which have been referred to as IntraDay Variability 
(IDV). It is known that IDV affects 20% to 50% of th e flat- 
spectrum radio source population (see Quirrenbach et al. 119921 
Lovell et al. 2008). The origin of the variability, however, is 
not completely understood. Assuming that the variability is in- 
trinsic to the source, causality arguments constrain the maxi- 
mum size of the emission region using the observed variability 
timescales. For IDV sources, this frequently implies very high 
brightness temperatures, exceeding by far the inverse-Compton 
limit (10 12 K; see Kellermann & Paulinv-Toth [T969b . 

Propagation effects, such as InterStellar Scintillation (ISS), 
were suggested as a likely alternative scenario. In this model, the 
variability is caused by the scintillation of very compact compo- 
nents (of the size of tens of juas) in the IDV source by a scat- 
tering screen placed between the source and the observer. The 
radiation can be either focused or defocused by the screen, caus- 
ing the measured flux density to be a function of the line of sight 
to the observer. Owing to the Earth's revolution around the Sun, 
the observer's line of sight towards the screen changes contin- 
uously, causing the detected rapid variations in the flux density 
measurements. Consequently, the faster the Earth's movement 
with respect to the scattering screen, the shorter the variability 
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timescale. Since the Earth's velocity follows a one-year periodic 
cycle, the relative velocity between the scattering screen and the 
observer should change accordingly, resulting in a seasonal cy- 
cle of the variability timescale. 

The detection of annual modulation patterns has been es- 
sential in demonstrating the source-extrinsic origin of IDV for 
the most extreme sources (timescales of hours, amplitudes of 
-100%), the so-called fast scintillators, namely PKS 0405-385, 
J1819 +384, and PKS 1257-326 (se e Ked ziora-Chudcz er et al . 
[19971 Dennett-Thorpe & de Bruvn 12000] Bignall et al. |2003). 
For type-II IDV sources (characterized by less extreme variabil- 
ity, with amplitudes of ~ 1 — 20% and timescales < 2 days, see, 
e.g., Quirrenbach et al. 1992), the situation is less clear. It may 
partially depend on their relatively long variability timescales, 
which make seasonal cycles harder to identify, given the unre- 
alistically long observing time that would be required to mea- 
sure accurate timescale estimates. Seasonal cycles have been de- 
tected in J 1128+592 (see Gabanyi et al. 120071 ) and claimed for 
other sources, such as S4 0917+62 (see Rickett et al. 2001 ). The 
MASIV survey (Lovell et al. l2003l and l2008] | showed that a sig- 
nificant part of the flux density variations observed in compact 
radio sources are due to ISS. On the other hand, the discovery 
of correlated optical-radio IDV in the flux density of the type-II 
source S5 0716+71 suggested that its origin may be intrinsic to 
the source (see Quirrenbach et al. 11991) . 

The detailed analysis of the variability characteristics of IDV 
sources plays an essential role in understanding the ISS contri- 
bution to their variability pattern. From 2005 to the end of 2009, 
the Max-Planck-Institut fur Radioastronomie and the Urumqi 
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Observatory carried out a program of intensive monitoring of 
type-II IDV sources with the Nanshan 25m radio telescope, lo- 
cated about 70 km south of Urumqi city (China). The main aims 
of the project are: (i) monitoring the long-term IDV charac- 
teristics of a larger sample of IDV sources, looking for sea- 
sonal cycles in their variability patterns, and (ii) detecting new 
IDV sources. The main target of the monitoring program were 
the type-II sources AO 0235+164, S5 0716+71, S4 0917+62, 
S4 0954+65, Jl 128+592, and HB89 1 156+295. During the ob- 
servations, about sixty additional sources were sporadically ob- 
served while attempting to find new IDV candidates. 

In the present paper, we focus on our results obtained for 
S4 0954+65, which is a well-studied BLLac object at redshift 
Z = 0.368 whose IDV features have been investigated since the 
early '90s (Wagner et al. 11993) . The origin of the variability 
is still controversial — S4 0954+65 is the second source for 
which a radio-optical correlation has been claimed (Wagner et 
al. [1990). However, the detection of extreme scattering events 
(Fiedler et al. [19871 [19941 Cimo et al. 120021 and the possi- 
ble presence of a seasonal cycle in the variability timescales of 
the source, which was revealed by a previous IDV monitoring 
project (see Fuhrmann 2004 ), seem to indicate that at least part 
of the variability is caused by extrinsic mechanisms, and that 
more than one scattering screen may be interposed between the 
source and the Earth. 

In the next section, we present an overview of the observa- 
tions and the data calibration procedure. The analysis methods 
used to extract the variability characteristics of the light curves 
are presented in Sec. 3. The main results are shown in Sec. 4. In 
Sec. 5 and 6, we discuss our results and present a summary of 
the paper, respectively. 

2. Observation and data calibration 

The Nanshan 25m telescope is equipped with a 4.8 GHz single- 
horn dual-polarization receiver provided, along with the new 
telescope driving program, by the MPIfR (for more details, see 
Sun et al. I2006I The receiver's bandwidth is 600 MHz. The 
flux density measurements were performed in cross-scan mode. 
Each scan consists of eight sub-scans in perpendicular directions 
across the source position in order to estimate the pointing off- 
sets in the two scanning directions and correct for them. The raw 
data were processed at the MPIfR via TOOLBOX, a Python- 
based software package that subtracts baseline drifts and fits a 
Gaussian profile to each cross-scan. The amplitude of the pro- 
file provides an estimate of the source's flux density expressed 
in units of antenna temperature (K). 

The data calibration was performed in semi-automatic mode 
(for more details, see Marchili 120091 Marchili et al. 2010) . The 
quality of the individual sub-scans was checked by considering 
their pointing offset, the full-width at half-maximum (FWHM) 
of their Gaussian profile, and the error in their flux density 
measurement. Sub-scans that showed pointing offsets of more 
than 100", deviations from the average FWHM (i.e. ~ 600") 
larger than 150", or uncertainties in the Gaussian fit estimate 
larger than 10% of the amplitude were automatically discarded. 
Subsequently, a weighted average of the remaining sub-scans 
was obtained with a weight proportional to the inverse of the 
squared errors, separately for the two scanning directions . After 
the pointing-error correction, we averaged the mean sub-scan 
amplitudes in the two scanning directions, providing a single 
flux-density measurement per scan. 

In each observing session, about 50% of the time was ded- 
icated to measurements of primary and secondary calibrators, 
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Fig. 1. The variability curves of S4 0954+65 (black dots) in 
April 2006 (upper panel) and April 2008 (lower panel). The 
turquoise squares show the residual variability in the fluxes of 
the secondary calibrators S5 0836+71 and S4 0951+69. 



i.e., compact radio-sources with steep spectra and no evidence of 
IDV. The 5 GHz flux densities of the so-called primary calibra- 
tors are fairly constant in time and well-known. The secondary 
calibrators, instead, show variability on timescales of months. 
We used the primary and secondary calibrators to estimate the 
variations in the gain of the antenna due to changes in both its 
parabolic shape and the position of the focus. We used the sec- 
ondary calibrators — which were selected to be at small angu- 
lar separations from the target sources — to correct for possi- 
ble changes in the collected flux densities caused by instabili- 
ties in the antenna-receiver system or changing weather condi- 
tions. The data calibration was completed by the conversion of 
measured antenna temperatures (in K) to flux densities (in Jy). 
The time-interpolated conversion factors were derived by means 
of the known flux densities of the primary calibrators, namely 
NGC 7027, 3C 286, and 3C 48. The complete data-calibration 
procedure guarantees a high level of accuracy, which, based on 
the residual variability in the calibrators, can be estimated in the 
range between 0.2 and 0.8%, under normal weather conditions. 

During the monitoring campaign, between August 2005 and 
December 2009, 41 IDV datasets for S4 0954+65 were col- 
lected. Two examples of light curves are shown in Fig.Q] In three 
cases, the observing sessions were too short (< 1.7d) to provide 
a reasonable estimation of the variability characteristics of the 
sources, were consequently discarded from the analysis we de- 
scribe in the following. An overview of the main characteristics 
of the remaining 38 sessions is reported in Table Q] In the few 
cases in which the duration of the observations is longer than 5 
days or shorter than 2.5 days, a bias of the estimated timescales 
toward, respectively, the slowest and the fastest variability com- 
ponents cannot be excluded. 



3. Variability analysis 

Following the scheme proposed by Quirrenbach et al. (2000 1 and 
Kraus et al. (12003ft . two different quantities were used to evaluate 
the significance and amplitude of the variability, namely the re- 
duced chi-squared test statistic x? an d the variability amplitude 
Y. 
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Table 1. Basic information about the 38 observing sessions at the 
Urumqi Observatory during which S4 0954+65 was observed. 
In Col. 1, we report the starting and ending dates of the experi- 
ments (format: dd.mm-dd.mm.yyyy),in Col. 2 the duration, and 
in Col. 3 the average number of measurements per hour for IDV 
sources (duty cycle). 
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The reduced chi-square (see, e.g., Bevington & Robinson 
1992 ) was used to test for given datasets the likelihood of the 
hypothesis of a constant flux density 



Xr 



= _J_ A / S;-<S > \ 2 



(1) 



where N is the number of datapoints, < S > the average flux 
density, 5; the i-th flux density measurement, and 6S\ its er- 
ror. A dataset was regarded as variable if the Xr 2 test returned 
a probability of constant flux below 0.1%, which is equivalent to 
a 99.9% probability for significant variability. 

An estimation of the variability in a dataset was provided by 
the modulation index, = 100 cr s /< S >, where cr s is the 

root mean square (rms) in flux density. However, this parameter 
does not take into account the residual noise in the calibrated 
data, which varies from epoch to epoch, mainly depending on 
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Fig. 2. Upper panels: the structure function plots of the April 
2006 (left panel) and April 2008 (right panel) datasets; the blue 
vertical lines show the estimated timescales. Lower panels: the 
wavelet-based analysis for the same epochs; the green vertical 
lines show the estimated timescales. The existence of two bumps 
in the right panels is the signature of two distinct variability com- 
ponents in the dataset. 



weather conditions. We defined the variability amplitude Y as 
an equivalent to a 3-cr confidence interval for Gaussian noise, 
modified by the quadratic subtraction of the rms noise floor, mo, 
to take into account the noise bias 



F = 3 



fit 



(2) 



3. 1 . Variability timescale 



The characteristic timescales of the variable datasets were esti- 
mated using two different methods of time series analysis — a 
first-order structure function (see Simonetti et al. II 9851 and a 
wavelet-based analysis. 

The structure function of a given time series {S (f,)},- is calcu- 
lated as 



SF(t) = -Y[S(.td-S(tj)f 



(3) 



where the sum is extended to the n pairs (f,-, tj) for which 
|| tj - tj ||=: r. Ideally, when a structure function is applied 
to a red-noise-like signal, it shows a monotonic increase that can 
be described as a power law. At the time lag r S f correspond- 
ing to the maximum coherency time of the signal, the struc- 
ture function shows a plateau, for which r s f provides an esti- 
mate of the characteristic timescale of the signal. Owing to both 
the uneven sampling and the noise in the analyzed data, the po- 
sition of the plateau is subject to some degree of uncertainty. 
For each plateau, it is possible to determine a minimum and a 
maximum timescale — respectively T m j n and r max . The char- 
acteristic timescale of the structure function r s f was calculated 
to be (r max + T m ; n )/2. The uncertainty in the timescale estima- 
tion was evaluated as the maximum between two terms, namely 
(i"max - T m ; n )/2 (which reflects the significance of the plateau) 
and 0.35 • (r sf 3/2 )/(ofcs d ) 1/2 (see Marchili et al. I20TTV where 
obsd indicates the duration of the observation. The latter term ac- 
counts for the significance of the timescale, which is limited by 
the small number of cycles that it is possible to observe, given the 
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limited duration of the observations. To estimate the proportion- 
ality factor, we calculated all the time-intervals between a local 
minimum (or maximum) and the following maximum (or mini- 
mum) for a set of light curves characterized by rapid variability; 
for each light curve, we calculated the average of the observed 
time-intervals (which provides an estimate of the timescale that 
is consistent with r S f) and their standard deviation (which pro- 
vides an estimate of its uncertainty). A linear regression of the 
uncertainties, plotted against the values of (r s t 3 ^ 2 ) /(obsd) 1 for 
the corresponding light curves, returned the proportionality fac- 
tor 0.35. 

The structure function in some cases was found to have 
more than one plateau, which is indicative of multiple variabil- 
ity timescales. In these cases, we identified the characteristic 
timescale with the shortest one. 

The wavelet-based algorithm that we developed is based on 
the Ricker ('Mexican hat') mother wavelet 



V(t) = c(l-t 2 )exp(- t -) 



(4) 



where C is a normalization factor and t is the time parameter. 
An extensive description of this analysis method can be found in 
Appendix A. The timescales obtained by applying the wavelet 
analysis are indicated with r wt . 

Examples of structure functions and wavelet results are pro- 
vided in Fig. |2]for the light curves shown in Fig. Q] 

4. Results 

The variability characteristics of S4 0954+65 are summarized in 
Table[2] Our study had two main outcomes, which are discussed 
below: (i) we found strong evidence of an abrupt change in the 
variability characteristics of S4 0954+65 between February and 
March 2008; (ii) the IDV timescales of the source seem to un- 
dergo an annual modulation. 

4.1. Change in the overall variability characteristics 

In the long-term light curve of the source (Fig. [3] upper panel), 
two considerably different phases can be recognized. From 
August 2005 up to February 2008, the flux density variations 
were generally moderate in intensity (the average peak-to-peak 
variation is ~ 0.2 Jy) and occurred on timescales of three 
months; March 2008 marked the starting point of a strong out- 
burst phase, which up to December 2009 was still ongoing, char- 
acterized by large variations (~ 0.4 Jy) on a timescale of about 
seven months. For simplicity, we refer to the first and the sec- 
ond timespan, respectively, as the pre-outburst and the outburst 
phase of S4 0954+65. 

The change in the long-term activity of the source was also 
reflected in the remarkable variation in its IDV characteristics. 
There are three lines of evidence for such a variation: 

- During the pre-outburst phase of the source, the IDV 
timescales provided by the wavelet-based method and the 
structure function (see Fig. [3] middle and lower panels) were 
highly correlated with each other. A linear regression applied 
to the timescales (Fig. |4] left panel) results in a correlation 
coefficient of 0.82 and a slope of 1.05, which proves the 
consistency in the definition of timescales used for the two 
methods. During the outburst phase, a similar analysis pro- 
vides a correlation coefficient of -0.11 with a slope of -0.19 
(Fig. H] right panel), implying that there is a clear lack of 



Table 2. The variability characteristics of S4 0954+65 in 38 
epochs of observation with the Urumqi Telescope. For each ses- 
sion (Col. 1), we report the average flux density (Col. 2), the 
variability amplitude (Col. 3), the reduced chi-square (Col. 4), 
the wavelet and structure function timescales (Col. 5 and 6, re- 
spectively), and the number of data points in each light curve 
(Col. 7). 
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1.01 


3.39 


3.05 


1.0+0.2 


0.4+0.2 


43 


2009.08.21 


1.07 


6.00 


9.57 


0.9+0.1 


1.3+0.3 


87 


2009.09.22 


1.22 


6.10 


11.02 


0.9+0.1 


1.5+0.4 


121 


2009.10.09 


1.22 


1.60 


1.38 






56 


2009.11.22 


1.27 


4.97 


5.38 


0.8+0.1 


1.0+0.2 


65 


2009.12.11 


1.27 


4.44 


6.51 


1.7+0.4 


0.5+0.2 


113 



correlation. It is useful to apply a linear regression to a com- 
parison sample comprising all the timescales estimated for 
the other main sources in the Urumqi monitoring program 
(i.e., AO 0235+164, S5 0716+71, S4 0917+62, Jl 128+592, 
and HB89 1156+295). This returns a correlation coefficient 
of 0.86 with a slope of 1.01, which proves the high de- 
gree of compatibility between the structure function and the 
wavelet-based method results. This demonstrates that the 
variability characteristics of the source are atypical during 
the outburst phase, as illustrated by a Kolmogorov-Smirnov 
test for the comparison of datasets. There is a 74% proba- 
bility that the distribution of the values of r s f/r wt during the 
pre-outburst phase is equal to the one estimated for the com- 
parison sample. A similar test for the r s f/r wt values during 
the outburst phase returns a probability below 0.1%. 
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Fig. 3. The long-term variability of S4 0954+65 (upper panel). In 
the middle panel, the wavelet timescales (red dots) and a three- 
point running average (blue line), showing the sudden increase 
starting from March 2008. The structure function timescales are 
shown in the lower panel. The red vertical line divides the ob- 
serving sessions into two phases characterized by remarkably 
different variability characteristics. 



- The wavelet timescales during the pre-outburst phase of 
S4 0954+65 follow a slowly decreasing trend, pinpointed by 
the three-point running average in the middle panel of Fig. [3] 
(blue line). A sudden increase (from ~ 0.5 to ~ 1.5 d) occurs 
in March 2008, after which the estimated timescales remain 
considerably longer than before. It is again helpful to apply 
a Kolmogorov-Smirnov test when comparing datasets; the 
probability that the two distributions of timescales are equal 
is below 1%. 

- We calculated an average structure function from the light 
curves obtained during the pre-outburst phase, and compared 
it with the average one for the outburst phaseQ. The differ- 
ence between the two can be seen in Fig. [5] In the struc- 
ture function of the pre-outburst phase, there is considerably 
more power on short timescales than in the structure func- 
tion of the outburst phase. On longer timescales, we find the 
opposite situation, suggesting that there is some source of 
slower variability that is more dominant during the outburst 
phase. 

Possible explanations of the large changes that occur in the 
IDV characteristics of S4 0954+65 at the beginning of its out- 
burst phase, are discussed in section [5] 

4.2. Annual modulation of the timescales 

We investigated the possible annual cycle in the variability 
timescales of S4 0954+65 by analyzing the evolution of the 
timescale r versus day of the year (DoY), separately for the 



1 Note that, in the case of ISS-induced variability, the rms of the flux 
density is expected to increase linearly with the average flux density 
(see Narayan 19921. This implies that the structure functions calculated 
during the outburst phase are expected to be systematically higher than 
the ones obtained during the pre-outburst phase. To make the results 
comparable, we applied the structure function to the light curves nor- 
malized by their average flux density. 



Slope: 1.01 (Corr. coeff.:0.86) 
Slope: 1.05 (Corr. coeff.: 0.82) 




Slope: -0.13 (Corr. coeff: -0.21) | 




[d] 



Fig. 4. Comparison between the structure function and 
wavelet results for the comparison sample (black dots) and 
the S4 0954+65 light curves. In the left panel, the source's 
timescales estimated up to February 2008 (cyan squares); in the 
right panel, the timescales from March 2008 to December 2009 
(green squares). The black, cyan, and green lines show the re- 
sults of applying linear regressions to the respective datasets. 
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Fig. 5. The average structure function calculated from the nor- 
malized light curves obtained during the pre-outburst phase 
(turquoise diamonds) and the outburst phase (red squares). 

wavelet-based and the structure function results. Following the 
scheme proposed by Qian & Zhang (2001), updated to the case 
of anisotropic scattering (see, e.g., Bignall et al. 2006. Gabanyi 
et al. l2007l l. we expressed r in terms of the vector s — which de- 
fines the orientation of the elliptical scintillation pattern — , the 
relative velocity between the scattering screen and the observer, 
v, the distance to the screen, D, and the anisotropy factor, r 



t(DoY) oc 



D- V? 



Vv 2 (DoY) + (r 2 - 1) (v(DoY) x s) 2 



(5) 



We note that the anisotropy may either be induced by an 
anisotropic scattering medium or be intrinsic to the source 
(caused by an anisotropic emitting component). 

We developed a code, based on a least squares fitting algo- 
rithm, to calculate the scintillation pattern parameters that most 
closely describe the detected variability timescales (for details, 
see Marchili 2009). We note that the distance cannot be esti- 
mated without making 'ad hoc' assumptions about the effective 
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angular size of the source; in the following, we adopt the model 
proposed by Qian & Zhang (2001), which assumes an intrinsic 
angular size of the scintillating component of 70 //as. In Section 
5.2, we modify this size estimate, to obtain an improved fit to the 
detected variation in the IDV pattern. 

The plots of the timescales versus DoY — the so-called 
annual modulation plots — , comprising all the estimated 
timescales between August 2005 and December 2009, are shown 
in Fig. |6l separately for the results provided by the two analy- 
sis methods. The wavelet results (left panels in the figure) sug- 
gest that two slow-down phases occur during the year, centered 
on May (DoY 140) and November (DoY 312), corresponding 
to a cycle of approximately six months. The structure function 
timescales (right panels in the figure) show a less clear cycle, 
with only two data-points being considerably higher than the 
others, corresponding to the observing sessions of December 
2005 and November 2006. Nevertheless, they appear to confirm 
that there is an annual cycle in the timescales that has a promi- 
nent slow-down phase peaking in November (DoY 320), and a 
much weaker one peaking in May (DoY 136) — in fair agree- 
ment with the wavelet results. The best-fit scintillation pattern 
parameters are reported in Table|3] The observed annual modula- 
tion cycle is compatible with the presence of a scattering screen 
at a distance between 180pc and 270 pc. The scintillation pattern 
is slightly anisotropic, with an anisotropy degree around 2 and 
an anisotropy angle around 90 deg. Were we to assume that the 
scattering screen is associated with the Ursa Major cloud com- 
plex, the estimated distance would be compatible with the value 
of 240 pc adopted by Heithausen et al. ( 119931) . but less consis- 
tent with the distance of only 100-120pc provided by Penprase 
(119931 ) for the molecular clouds MBM 29-MBM 31 in the com- 
plex. A distance to the screen of between 120 and 290 pc can 
also be deduced from the Hipparcos catalog (Perryman & ESA, 
1997; van Leeuwen. 120071 ). by looking at the parallaxes of the 
stellar objects within a two-degree circle around S4 0954+658. 



5. Discussion 

The existence of a break in the variability characteristics of the 
source, described in Section|4] raises the question of whether this 
break affects the search for an annual variation in the character- 
istic timescales. The consequences of the March 2008 break can 
be evaluated by separately analyzing the variability timescales 
before and after it occurs. The corresponding annual modula- 
tion plots and scintillation pattern parameters are presented in 
Fig. [7] and Table [3] for the wavelet and structure functions sepa- 
rately. In the upper panels of Fig. [7] the August 2005 - February 
2008 timescales are shown as magenta squares, while the March 
2008 - December 2009 timescales are shown as black circles. 
In the middle and lower panels, we show 40-day averages of the 
timescales during the pre-outburst and the outburst phases, re- 
spectively. In all plots, the green and turquoise lines display the 
best annual modulation fits to the data. Apparently, the yearly 
cycles found before and after March 2008 in the wavelet-based 
results are quite similar, with slow-down phases at DoY 144 and 
320 for the former and at DoY 137 and 309 for the latter. In the 
case of the structure function results, instead, the difference in 
the cycles detected in the two timespans is significant. To under- 
stand the nature of the change in the variability characteristics of 
S4 0954+65 between the pre-outburst and the outburst phases, it 
is essential to find the reason why the wavelet and the structure 
function results are affected differently. 



5. 1. Multiple timescales in the light curves after February 
2008 

The general increase in the wavelet timescales after February 
2008 could be explained in terms of a general slow-down in the 
variability of the source. This, however, could not account for 
the systematic discrepancies between the results of the two anal- 
ysis methods, highlighted in Fig.|4] These discrepancies are most 
likely the consequence of the appearance, starting from February 
2008, of multiple variability contributions in the light curves, of 
different timescales and comparable strengths. The wavelet and 
structure function plots in the right panels of Fig. [2] which show 
the results for one of the light curves collected during the out- 
burst phase, confirm this interpretation. The wavelet plot shows 
two bumps, corresponding to timescales of 2.2 and 0.5 days, 
which indicate that there are two distinct variability contribu- 
tions; the two plateaus in the structure function plot, which have 
similar timescales, can be explained in the same way. 

The discrepancies between the estimated timescales reveal 
the differences between the two analysis methods. The wavelet 
timescales are generally longer than the structure function ones, 
suggesting that the slowest variability contributions are the most 
powerful. For the structure function, it was our choice to iden- 
tify the characteristic timescale with the one corresponding to 
the first significant plateau of the function. This rigid approach, 
which was forced by the necessity to minimize the arbitrariness 
of the definition of the characteristic timescales, should not in- 
terfere with the detection of possible annual modulation cycles. 
Assuming that a fast and a slow timescales are both detected 
in the majority of the light curves — which is the case for the 
observing sessions during the outburst phase — and that they 
are caused by the same scattering screen, they should follow the 
same annual modulation cycle. Therefore, theoretically, this can 
be identified by systematically looking at the fastest variabil- 
ity component — which is the one that can be estimated with 
the highest accuracy, given the increase in the uncertainty with 
lengthening timescale. In practice, however, the situation is more 
complicated; the large difference in the annual modulation plots 
of the structure function timescales before and after March 2008 
demonstrates that a proper interpretation of the results cannot be 
achieved without simultaneously taking into account the infor- 
mation provided by both of the time analysis methods. 

5.2. The nature of the different variability components 

We present below two possible models for explaining the annual 
modulation plots of the outburst phase of S4 0954+658 (see Fig. 
[7] lower panels). Both assume that a source of variability is the 
propagation effects caused by a scattering screen between the 
source and the observer, and that a reasonable estimate of its 
main parameters is provided by the pre-outburst wavelet results 
(see Tab. [3] and Fig. [7] middle-left panel). These assumptions 
are based on the existence of an annual cycle in the pre-outburst 
timescales and the detection of a similar cycle in the wavelet 
results during the outburst phase of the source. 

In the first model, we hypothesize that the second contribu- 
tion to the variability is provided by a source-intrinsic mech- 
anism; the corresponding variability timescale should show no 
annual modulation effect. The frequent detection of timescales 
of about one day in the wavelet and structure function results 
suggests that the source-intrinsic process could be described in 
terms of a variability contribution of constant timescale. The 
present scenario is represented in Fig. [8] left panel. The wavelet 
and structure function results are plotted along with a one- day 
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Table 3. The scintillation pattern parameters deduced by the distribution of the timescales as a function of the day of the year. In 
Col. 1, we report the timespan taken into account for the estimation of the parameters; in Col. 2, the method used for estimating the 
timescales; in Col. 3 and 4, the relative velocity between the scattering screen and the observer, projected into right ascension and 
declination; in Col. 5, the distance to the screen; in Col. 6 and 7, the anisotropy degree and angle. 



Epochs 


Method 


V ra 


Vdec 


D 




r 


e 






[km/s] 


[km/s] 


[kpc] 






[deg] 


08.2005 - 12.2009 


WT 


10 ±3 


10 ±5 


0.27 ± 0.03 


2.3 


±0.3 


84 ±6 


08.2005 - 02.2008 


WT 


11 ±5 


8±7 


0.25 ± 0.04 


3.0 


±0.5 


90 ±6 


03.2008 - 12.2009 


WT 


11 ±5 


11 ±5 


0.27 ± 0.04 


1.7 


±0.3 


82 ± 12 


08.2005 - 12.2009 


SF 


12 + 6 


-1±5 


0.18 ±0.03 


1.8 


±0.3 


90 ± 10 


08.2005 - 02.2008 


SF 


12 ± 5 


-5 ±6 


0.23 ± 0.04 


1.7 


±0.3 


82 ± 12 


03.2008 - 12.2009 


SF 


5±9 


9±9 


0.14 ±0.08 


1.3 


±0.3 


33 ± 15 



constant value (green line) and an annual modulation pattern 
(blue line), which assumes an intrinsic angular size of the scin- 
tillating component of 50 yuas. To improve the fit to the data, 
the anisotropy degree of the scattering pattern was increased 
from three to seven. The most reasonable way of explaining 
this change is to hypothesize that the origin of the anisotropy is 
the emitting component of the source; structural changes in this 
component would also explain the observed break in the vari- 
ability characteristics of S4 0954+658. 

In the second model, both the variability contributions are 
caused by ISS through the same scattering screen, parameterized 
according to the fit to the pre-outburst wavelet results. The corre- 
sponding annual modulation curves, plotted in the right panel of 
Fig. [8] assume intrinsic angular sizes of the scintillating com- 
ponents of 90 (blue line) and 30 yuas (green line). Following 
Marscher et al. (11 9791 . and considering an amplitude of the IDV 
of ~ 0.1 Jy, from the size of the smaller component we can cal- 
culate its brightness temperature. This turns out to be higher than 
10 12 K; in order to avoid a violation of the inverse-Compton 
limit, we have to assume that the emitting component is mov- 
ing relativistically with a Doppler boosting factor > 3. 

The characteristic timescales estimated during the outburst 
phase of S4 0954+658 seem compatible with both models. The 
source-intrinsic scenario provides a simple explanation of the 1 1 
occurrences of timescales of about one day detected in the 17 
observing sessions, but its fit to the wavelet timescales is worse 
than the one obtained through the fully source-extrinsic scenario. 
None of the models explains the long timescale detected by the 
wavelet method in July 2008. 

The limited amount of information provided by the variabil- 
ity analysis does not allow us to reach any final conclusion. 
Whatever the origin of the additional variability component, its 
appearance is in all cases remarkable, especially when associ- 
ated with the outburst phase of the source. It seems unlikely that 
the simultaneity of the two events is accidental. Since the out- 
burst phases of blazars are often explained by changes in the 
structure of the sources, such as the emission of new emitting 
components (see, e.g., Mutel et al. [1990), and since the size 
and luminosity of the emitting regions are supposed to play 
a fundamental role in IDV, the existence of a correlation be- 
tween outburst phases and IDV characteristics of a source ap- 
pears pretty reasonable. A possible connection between changes 
in the variability patterns of IDV sources and variations in their 
structures was reported by Fuhrmann et al. (2002 1 and Macquart 
& de Bruyn (2007), and our findings seem to confirm this pic- 
ture. We note that many AGNs show intrinsic structural vari- 
ability on timescales of months to years. These varying source 
structures complicate the search for annual modulation patterns, 



which may be detectable in AGNs only during sufficiently long 
periods of quiescence. 

6. Summary 

We have reported on the evolution of the intraday variability 
characteristics of S4 0954+65 during the period from August 
2005 to December 2009. For this source, 41 IDV observ- 
ing sessions have been collected at the Urumqi Astronomical 
Observatory at a frequency of 4.8 GHz; 38 of these sessions 
were long enough to provide a reliable estimation of the vari- 
ability characteristics. 

The data analysis has been performed using standard tools of 
statistics and time series analysis, such as the variability ampli- 
tude, the chi-square test for evaluating the probability of a con- 
stant flux density and the first-order structure function. A new 
method for the estimation of the characteristic timescale — a 
modified wavelet analysis based on the Ricker mother wavelet 
— has been developed and successfully applied to the data. The 
results provided by this new algorithm show a high degree of 
correlation with the ones obtained using a structure function 
analysis. We collated a comparison data set of more than 120 
light curves, comprising the light curves of all the main targets 
of the monitoring program except S4 0954+65; for this data 
set, a linear regression between wavelet and structure function 
timescales recovers a correlation coefficient of 0.86 and a slope 
of 1.01, which demonstrates the consistency between the defini- 
tions of timescales in the two methods. The wavelet estimation 
of the timescales offers the important advantages of providing 
an unambiguous definition of the variability timescale and the 
possibility of full automation. 

The combined use of a wavelet and structure function anal- 
ysis proved to be a powerful tool for tracing the evolution of 
the IDV timescales of S4 0954+65. Both methods reveal that 
there is an annual modulation cycle in the variability timescales, 
suggesting an essential contribution of propagation effects to the 
detected variability. The seasonal cycle is not as strong as the 
ones detected in fast scintillators, such as J1819+3845 (Dennett- 
Thorpe & de B ruyn |2000T >. or in the type-II source Jl 128+5925 
(Gabanyi et al. 2007), but the agreement between the cycles ob- 
served in the wavelet timescales before and after February 2008, 
with slow-down phases around DoY 140 and 310, and the con- 
firmation of a slow-down phase around DoY 320 in the struc- 
ture function results are strong arguments in favor of a source- 
extrinsic origin of the variability. 

The comparison of the results from the two methods con- 
tributed to the discovery of a substantial change in the IDV char- 
acteristics of the source, occurring between February and March 
2008 at the beginning of a strong outburst phase. The signature 
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Fig. 6. Annual modulation plots based upon the wavelet results 
(left panels), and the structure function results (right panels). The 
upper panels show all the estimated timescales, while lower pan- 
els show monthly averages. The green and turquoise lines show 
the annual modulation patterns which best fit the data. 

4- 




180 
DoY 

Fig. 7. Annual modulation plots, as in Fig.[6j in the upper pan- 
els, the timescales obtained for the epochs between August 2005 
and February 2008 are highlighted in magenta squares. In the 
middle and lower panels are displayed the annual modulation 
plots (green and turquoise lines) and the 40-day timescales aver- 
ages for the August 2005 - February 2008 and the March 2008 - 
December 2009 epochs. The red vertical lines indicate the peaks 
in the annual modulation patterns obtained from the complete 
set of timescales. 



of this change is a sudden lengthening in the average timescales 
detected by the wavelet-based analysis, an increase in the vari- 
ability amplitude measured through the structure function, and 
a dramatic weakening in the correlation among the results pro- 
vided by the two methods. The variability detected in the light 
curves is likely the result of the superposition of two independent 
contributions. The data seem to favor a source-extrinsic origin 
for both the contributions, but the results cannot be unambigu- 
ously interpreted. Nevertheless, the clear correlation between the 
long-term flux density variations in S4 0954+65 — likely asso- 




180 
DoY 



360 



180 
DoY 



Fig. 8. The wavelet (black dots) and structure function (red 
squares) timescales calculated during the outburst phase of 
S4 0954+658 are indicative of multiple variability contributions. 
In the left panel, the variability is interpreted as the superposi- 
tion of a source-extrinsic effect (blue line) and a source-intrinsic 
process (green line). In the right panel, both the variability con- 
tributions are attributed to ISS and the two annual modulation 
patterns have the same scattering screen parameters, the only dif- 
ference being the size of the scintillating components (90 yuas for 
the blue line, 30 //as for the green line). 



ciated with changes in the source structure — and its IDV fea- 
tures shows that the intrinsic properties of the source exert a sig- 
nificant influence on the characteristics of the variability on the 
shortest timescales. 
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Appendix A: Wavelet analysis 

The wavelet-based analysis we used to estimate the timescales 
in our datasets is based on the Ricker mother wavelet 

V(t) = c(l-t 2 )exp(-^), (A.l) 

where C is a normalization factor and t is the time. Given the 
scale s and the time-shift /, we define the wavelet 

™-£(>-^M-^) 

i 

and the wavelet transform 

«s, = f (St- < S » (l - <*Z£) exp( - (A3) 

i 

For each variability curve, we calculated if/(s, I) on a range of 
scales s, from twice the average sampling to half the total length 
of the curve, with / varying between (to + s/2) and (f# - s/2); 
we defined the variability timescale as V3 %, where sm is the 
scale for which I) is maximum, noting that, given a scale s, 
V3 s is the time lag between the minimum and the maximum of 

The results provided by the method described above are not 
satisfactory. To allow a proper comparison between the wavelet 
transforms, Yjf ^s, i(td should be close to zero for each set of 
parameters (s, I). This condition is often unfulfilled, because of 
the uneven sampling of the variability curves and their limited 
duration. To correct for this, we introduced a modified wavelet 
transform 



As for the standard wavelet analysis, the variability timescale 
is defined as V3 sm, where sm is the scale that corresponds to 
the maximum value of <p s j. The standard wavelet analysis is 
generally more sensitive to slow variability components than 
to fast ones, often detecting unreasonably long characteristic 
timescales. The results provided by the modified wavelet, in- 
stead, appear to be much more reliable. 

The modified wavelet transform has a strong advantage over 
the structure function analysis in dealing with light curves with 
variability on multiple timescales. While the identification of the 
characteristic timescale in the structure function requires a sub- 
jective interpretation of which timescale is the most predomi- 
nant, the wavelet analysis always provides a unique value, be- 
cause the set (s, I) that maximizes <f> s j is always univocally de- 
fined. This completely removes the ambiguity when interpreting 
the analysis results — regardless of who performs the analysis, 
the same timescale is found. It also implies that the wavelet anal- 
ysis can be fully automated, which can be very useful in dealing 
with large datasets. On the other hand, one should always keep in 
mind that the timescale obtained by the wavelet method, strictly 
speaking, does not provide an estimation of the characteristic 
variability timescale. It refers, instead, to the strongest variabil- 
ity feature in the light curve. When the variability timescale in a 
light curve is much shorter than the duration of the observations, 
the structure function should be regarded as a more suitable tool 
for the data analysis. 

We note that the intraday variability of blazars — such as 
S4 0954+65 — at radio frequencies can be generally described 
as red noise. The amplitude of the variability is stronger on 
longer timescales, which implies that the variability timescale is 
typically on the same order of magnitude as the duration of the 
observations. In these conditions, the identification of the char- 
acteristic timescales with the strongest variability feature in the 
light curves appears quite reasonable. This conclusion is con- 
vincingly supported by the high degree of correlation between 
the structure function and the wavelet results (see Section[4]i. 

The error in the timescale estimates provided by the modified 
wavelet transform can be ascribed to two main contributions. 
The first consists in the combined effect of the uneven sampling 
and the uncertainty in the flux density measurements; the second 
has to do with the limited duration of the observations, which 
causes a gradual decrease in the reliability of the wavelet trans- 
forms on larger scales s. Both effects are reflected in the width 
of the peak \jj(sM,l)- Wider peaks are indicative of a low con- 
trast among the strengths of the variability on different scales s, 
hence larger errors in the estimation of the timescale. We per- 
formed tests with synthetic light curves of different durations, 
sampled both evenly and unevenly. We found that a reasonable 
evaluation of the error in the timescales is given by (52 - *i)/2, 
where S2 and S\ are respectively the maximum and the minimum 
scales for which i[/(s, I) equals 0.95 • ^(sm, 0- 



4>(s,l) = C s jif/(s,l), 
where C s / is defined as 
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